import numpy as np
import matplotlib.pyplot as plt
import time
def f(x):
    return np.sqrt(2/np.pi)*np.exp(-x**2/2)
def MH_normal(N):
    d = np.zeros(N)
    x = 0
    for i in range(N):
        y = np.random.rand() * 2 - 1 + x
        h = np.min([1, f(y) / f(x)])
        U = np.random.rand()
        if U < h:
            x = y
        d[i] = x
    return d
N=100000
def MH_view(N):
    k, X = N,MH_normal(N)
    x = np.linspace(-6,6,10000)
    plt.hist(X[0:k], bins=50,density=1, label="Stat");
    plt.plot(x,1/np.sqrt(2*np.pi)*np.exp(-x**2/2), label=r"PDF")
    plt.legend()
    plt.savefig('MH.png')
MH_view(N)
X=np.arange(1000,100000,1000)
Y=[]
for x in X:
    t_0=time.time()
    MH_normal(x)
    t=time.time()-t_0
    Y.append(t)
plt.plot(X,Y,'r.')   
plt.xlabel(r"size of Normal Distribution Samples")
plt.ylabel(r"time")
plt.savefig('MH_time.png')